#ifndef _VOLUME_REGION_H_
#define _VOLUME_REGION_H_


class VolumeRegion{
public:
	virtual ~VolumeRegion(){};
	virtual bbox_t WorldBound()const=0;
	virtual bool IntersectP(const ray_t& r, float *t0, float *t1)const = 0;
	virtual spectrum_t sigma_a(const point3f_t &, const vector3f_t &)const = 0;
	virtual spectrum_t sigma_s(const point3f_t &, const vector3f_t &)const = 0;
	virtual spectrum_t lve(const point3f_t &, const vector3f_t &)const = 0;
	virtual float p(const point3f_t&, const vector3f_t &, const vector3f_t&)const = 0;
	virtual spectrum_t tau(const ray_t &ray, float step , float offset)const = 0;

	virtual spectrum_t sigma_t(const point3f_t &, const vector3f_t &wo)const;
};


float phase_hg(const vector3f_t &w, const vector3f_t &wp, float g);
float phase_schlick(const vector3f_t &w, const vector3f_t &wp, float g);
template<typename Derived>
class DensityRegion:public VolumeRegion{
public:
	DensityRegion(const spectrum_t &sa, const spectrum_t &ss, float g,
		const spectrum_t& emit, const transform_t &volume_to_world)
		:sig_a_(sa),sig_s_(ss), le_(emit),g_(g){
			transform_inverse(world_to_volume_,volume_to_world);
	}

	float Density(const point3f_t &pobj)const{return derived().Density(pobj);}

	spectrum_t sigma_a(const point3f_t &p, const vector3f_t &)const
	{
		point3f_t pobj;
		transform_point(pobj, world_to_volume_,p);
		spectrum_t ret;
		vsmul(ret, (Density(pobj)) , sig_a_);
		return ret;
	}
	spectrum_t sigma_s(const point3f_t &p, const vector3f_t &)const
	{
		point3f_t pobj;
		transform_point(pobj, world_to_volume_,p);
		spectrum_t ret;
		vsmul(ret,(Density(pobj)) , sig_s_);
		return ret;
	}
	spectrum_t sigma_t(const point3f_t &p, const vector3f_t &)const
	{
		point3f_t pobj;
		transform_point(pobj, world_to_volume_,p);

		spectrum_t ret;
		vadd(ret,sig_a_, sig_s_);
		vsmul(ret,(Density(pobj)) , ret);
		return ret;
	}
	spectrum_t lve(const point3f_t &p, const vector3f_t &)const 
	{
		point3f_t pobj;
		transform_point(pobj, world_to_volume_,p);
		float d = Density(pobj);
		spectrum_t ret;
		vsmul(ret,d,le_);
		return ret;
	}
	float p(const point3f_t& p, const vector3f_t &w, const vector3f_t &wp)const{
		return phase_hg(w,wp,g_);
	} 
	spectrum_t tau(const ray_t& r, float step_size, float offset)const;
private:
	Derived& derived(){return static_cast<Derived>(*this);}
	const Derived& derived()const{return static_cast<Derived>(*this);}
protected:
	spectrum_t sig_a_, sig_s_, le_;
	float g_;
	transform_t world_to_volume_;
};


template<typename Derived>
spectrum_t DensityRegion::tau(const ray_t& r, float step_size, float offset)const
{
	float t0,t1;
	float length = std::sqrt(vdot(r.d,r.d));
	spectrum_t tau;
	vclr(tau);
	if(length == 0.f)return tau;
	ray_t rn;
	rn = r;
	rn.o = r.o;
	rn.d.x = r.d.x/length;
	rn.d.y = r.d.y/length;
	rn.d.z = r.d.z/length;
	rn.mint = r.mint * length;
	rn.maxt = r.maxt * length;
	if(!IntersectP(rn, &t0, &t1))return tau;
	t0 += offset * step_size;
	spectrum_t vtmp;
	while(t0 < t1)
	{
		point3f_t p;
		vector3f_t d;
		rpos(p,rn,t0);
		vneg(d,rn.d);
		vtmp = sigma_t(p,d);
		vadd(tau,tau,vtmp);
		t0 += step_size;
	}
	vsmul(tau, step_size, tau);
	return tau;
}

float phase_hg(const vector3f_t &w, const vector3f_t &wp, float g)
{
	float cos_theta = vdot(w,wp);
	return 1.f/(4.f*PI) * (1.f - g * g)/powf(1.f + g*g-2.f*g*cos_theta, 1.5f);
}
float phase_schlick(const vector3f_t &w, const vector3f_t &wp, float g)
{
	float alpha = 1.5f;
	float k = alpha * g -(1.f-alpha)*g*g*g;
	float kcostheta = k * vdot(w,wp);
	return 1.f / (4.f * PI) * (1.f - k*k)/((1.f-kcostheta)*(1.f-kcostheta));
}


class VolumeGridDensity:public DensityRegion<VolumeGridDensity>
{
	typedef DensityRegion<VolumeGridDensity> base_type;
public:
	VolumeGridDensity(const spectrum_t &sa, const spectrum_t &ss, float g,
		const spectrum_t &emit, const bbox_t &bound, const transform_t& volume_to_world,
		int x, int y, int z, const float *dense)
		:base_type(sa,ss,g,emit,volume_to_world),nx_(x),ny_(y),nz_(z),bound_(bound)
	{
		density_ = new float[x*y*z];
		::memcpy(density_,dense,nx*ny*nz*sizeof(float));
	}
	~VolumeGridDensity(){delete []density_;}

	float Density(const point3f_t& pobj) const;
	bool IntersectP(const ray_t& r, float *t0, float *t1)const;
private:
	float *density_;
	const int nx_,ny_,nz_;
	const bbox_t bound_;
};
#endif